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Abstract 

This  paper  focuses  on  the  thermal  engineering  design  and  analysis  of  solid  oxide  fuel  cell  (SOFC)  units,  with  emphasis  on  cell  performance 
and  component  design.  In  engineering  practice,  insulation  materials  would  be  deployed  as  the  enclosure  of  an  SOFC  stack  to  reduce  the  heat  loss 
to  the  environment.  In  this  work,  a  computational  methodology  has  been  implemented  to  characterize  the  thermal  engineering  performance  of 
a  planar  SOFC.  The  present  calculation  procedure  integrates  the  steady-sate  electrochemical  reactions  of  the  SOFC  with  finite-element  models 
for  thermo-mechanical  analyses  of  the  interconnect  through  iteration  processes,  so  that  a  unified  temperature  distribution  with  heat  loss  effect 
can  be  obtained.  Present  results  show  that  the  convergent  rate  of  the  adopted  methodology  is  quite  efficient,  and  that  the  temperature  patterns 
are  compatible  with  those  reported  in  the  literature.  Furthermore,  this  work  has  also  developed  a  bulk  heat-transfer  model  for  simplified  design 
analysis.  The  concept  of  total  heat  resistance  is  employed  to  facilitate  the  one-dimensional  (ID)  analyses  and  to  determine  the  predominant 
parameters  that  affect  heat- transfer  behaviour.  Moreover,  some  accommodation  factors  have  been  deduced  to  correlate  the  ID  results  of  lateral 
heat  transfer  with  those  of  two-dimensional  (2D)  finite-element  analyses,  as  this  will  be  beneficial  for  rapid  prototyping  processes. 

©  2004  Elsevier  B.Y.  All  rights  reserved. 
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1.  Fundamentals  and  approach 

Among  the  various  types  of  fuel  cell,  the  solid  oxide  fuel 
cell  (SOFC)  belongs  to  the  high- temperature  category;  its  op¬ 
erating  temperature  ranges  from  about  600  to  1000  °C.  The 
working  principle  of  the  SOFC  is  illustrated  schematically  in 
Fig.  1.  Fuel  is  fed  into  the  anode  side,  while  on  the  opposite 
cathode  side,  air  is  typically  utilized  as  the  oxidant.  The  oxi¬ 
dant  generates  02~  ions  that  serve  as  charge  carriers  through 
the  electrolyte  to  the  anode,  where  they  react  with  fuel  to 
release  electrons.  By  guiding  electrons  via  an  outer  electric¬ 
conducting  loop  back  to  the  anode,  this  device  generates  di¬ 
rect  current  (dc)  electricity.  One  of  the  major  advantages  of 
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the  SOFC  is  its  fuel  flexibility.  For  example,  the  fuel  can 
be  a  mixture  consisting  of  hydrogen  (FL),  carbon  monoxide 
(CO),  methane  (CEL)  and  even  some  higher  hydrocarbons. 
This  feature  avoids  the  expense  of  producing  hydrogen  of 
high  purity  that  is  demanded  by  low-temperature  fuel  cells. 

The  basic  cell  consists  of  a  positive  electrode 
(cathode)-electrolyte-negative  electrode  (anode)  (PEN)  as¬ 
sembly  and  two  end-plates,  with  piping  on  each  side.  Flow 
channels  for  fuel  and  oxidant  are  integrated  on  the  end  plates 
at  the  anode  and  cathode  sides,  respectively.  In  practical  ap¬ 
plications,  multiple  cells  are  piled  up  to  form  a  stack  that 
corresponds  to  serial  connections  of  the  electric  loop  for  an 
individual  cell.  In  this  case,  interconnects  with  double-sided 
flow  channels  will  lay  between  the  PENs.  A  scheme  of  a  stack 
assembly  is  shown  in  Fig.  2. 

To  design  and  optimize  a  stack  configuration,  it  is 
necessary  to  establish  analytical  tools  for  simulating  the 
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Nomenclature 

Symbols 

A 

contact  area  in  channel  for  convection 

C 

molar  flow  rate 

Cp 

specific  heat 

F 

Faraday  constant 

G 

Gibbs  function 

Gr 

Grashof  number 

g 

gravity 

H 

height  of  channel 

h 

heat-transfer  coefficient 

I 

current 

k 

thermal  conductivity 

L 

characteristic  length 

Nu 

Nusselt  number 

P 

pressure 

Pr 

Prandtl  number 

Q 

heat  source 

q 

heat  flux  per  unit  surface 

R 

gas  constant,  heat  resistance 

Re 

Reynolds  number 

T 

temperature 

U 

total  heat-transfer  coefficient 

u 

velocity 

V 

voltage 

w 

width  of  channel 

X 

thickness 

Greek  letters 

p 

coefficient  of  expansion 

A 

difference 

V 

kinematic  viscosity 

X 

oxygen  to  water  ratio 

Subscripts  and  superscripts 

0 

standard  state 

a 

air 

act 

activation 

con 

concentration 

f 

fuel 

i 

species 

L 

characteristic  length 

m 

mean 

nst 

Nernst 

0 

outer 

ohm 

ohmic 

op 

operation  state 

pol 

polarization 

s 

solid 

v  direction 

y 

y  direction 

w 

wall 

FUEL  (CH.  or  FL  or  CO) 


FL  CH 


o2,  n2 


Electrolyte 


Cathode 


T 


AIR  (CL) 


Anode 

B2  +  O2  _ H20  +  2e- 

CO  +  O2  —  C02  +  2e- 


Cell  Overall 


H2  +  l/202-H20 
CO  +  l/202  — C02 


* 


+ 

Cathode 
02+4e-— 202 


Fig.  1.  Schematics  of  operating  principles  of  a  SOFC. 


cell  performance.  Reports  of  the  modelling  of  planar 
SOFCs  have  been  increasing  during  the  past  few  years, 
and  most  of  these  have  dealt  with  steady  operation.  Some 
previous  studies  will  be  briefly  reviewed  here  to  illustrate 
the  evolution  of  modelling  developments.  Pioneer  work  [1] 
examined  a  three-dimensional  (3D)  model  for  steady  and 
load-following  operations  to  analyze  the  cell  characteristics 
as  well  as  the  distribution  of  temperature,  current  density 
and  gas  species  under  the  influence  of  various  parameters. 
The  basic  assumptions  included  uniform  flow  in  the  cell, 
no  pressure  drop  and  adiabatic  conditions  on  the  periphery 
of  stack.  These  have  been  adopted  in  some  later  studies. 
Iwata  et  al.  [2]  investigated  the  effects  of  fuel  re-circulation 
ratio,  operating  pressure,  flow  configuration  and  thermal 
boundary  conditions  on  the  distribution  of  temperature  and 
current  density.  The  conservation  equations  were  written  in 
discretized  forms  for  cell  elements,  which  were  arranged 
periodically  in  span- wise  direction  to  complete  a  unit  cell. 
Similar  assumptions  as  those  in  [1]  were  implemented.  In 
addition,  each  segmented  cell  was  assumed  to  have  one  repre¬ 
sentative  temperature  and  an  additional  option  of  a  radiation 
boundary  could  be  set  at  the  upper  and  lower  surfaces. 

Some  investigators  have  employed  commercial  computa¬ 
tional  fluid  dynamics  (CFD)  codes  to  perform  thermofluid 
analyses  of  planar  cells,  in  which  the  computational  domain 
has  been  focused  on  either  flow  channels  [3,4]  or  cell  compo¬ 
nents  [5].  For  example,  Yakabe  et  al.  [3]  presented  the  distri¬ 
bution  of  chemical  species,  temperature,  potential  and  current 
density  in  flow  channels  of  a  planar  SOFC.  Their  model  was 
based  on  a  simplified  single-unit  cell  with  a  parallel  flow  pat¬ 
tern  and  a  finite  volume  that  simulated  half  of  the  repeating 
unit  located  in  the  middle  part  of  a  one-cell  stack.  The  ther¬ 
mal  boundaries  employed  adiabatic  and  cyclic  conditions  at 
the  edges  and  the  surfaces  connecting  to  the  adjacent  units  in 
a  span-wise  direction,  respectively.  For  the  sake  of  simplic¬ 
ity,  radiation  losses  were  neglected  in  the  steady-state  energy 
equation,  and  the  diffusion  overpotential  was  also  omitted 
in  the  discrete  electrochemical  model  (ECM).  Formulae  for 
radiant  heat-exchange  inside  the  channels  were  derived  in 
follow-up  work  but  were  implemented  only  for  one  calculat¬ 
ing  condition  because  of  the  amount  of  time  required,  which 
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Fig.  2.  Schematics  of  SOFC  stack  assembly. 


turned  out  to  be  10  times  that  of  the  non-radiant  counter¬ 
part.  The  effect  of  radiation  has  been  considered  to  be  small, 
however,  due  to  the  large  aspect  ratio  of  the  flow  channels 

[4] . 

An  alternative  approach  was  presented  by  Recknagle  et  al. 

[5] ,  in  which  3D  model  geometries,  including  internal  man¬ 
ifolds,  PENs,  interconnects  and  glass  seals,  were  created  to 
simulate  generic  stack  designs,  either  for  cross-  or  parallel- 


flow  configurations.  Individual  flow  channels  created  by  webs 
in  the  interconnect  were  not  modelled;  instead,  the  region  was 
assumed  to  be  one  large  flow  passage.  The  effect  of  radiation 
was  perceived  to  be  small  relative  to  the  other  modes  of  heat 
transfer  and  was  neglected.  Cyclic  boundary  conditions  were 
imposed  at  the  top  and  bottom  of  the  model  domains,  while 
the  lateral  walls  were  assumed  to  be  adiabatic.  The  above 
thermo-fluid  model  was  combined  with  a  unit-cell  electro¬ 
chemistry  model  to  form  a  simulation  tool  for  a  planar  SOFC. 

Due  to  the  high  operating  temperature  of  SOFCs,  heat  loss 
to  the  environment  becomes  inevitable  and,  in  turn,  will  affect 
the  overall  system  performance.  Heat  transfer  across  the  walls 
at  the  periphery  of  the  stack  is  hard  to  estimate,  however,  so 
that  this  type  of  thermal  boundary  has  been  ignored  in  the  en¬ 
ergy  equation  in  many  reports  of  the  characterization  of  cell 
performance,  including  electrochemical,  thermal  and  flow  as¬ 
pects.  As  mentioned  by  some  investigators,  it  is  difficult  to 
determine  the  exact  boundary  at  the  edge  parts.  Hence,  all  the 
aforementioned  studies  have  employed  adiabatic  boundaries 
for  the  energy  equation  on  the  periphery  of  the  computational 
domain,  regardless  of  flow  channels  or  cell  components. 

In  this  work,  an  efficient  computational  methodology  is 
proposed  to  deal  with  the  behaviour  discussed  above,  of 
which  the  procedures  and  some  preliminary  results  have  been 
presented  to  local  and  international  scientific  communities 
[6,7].  At  present,  a  computer  code  for  analyzing  the  steady- 
sate  electrochemical  reactions  has  been  developed,  by  which 
the  characteristics  of  SOFC  can  be  calculated.  On  the  other 
hand,  finite-element  models  for  thermo-mechanical  analy¬ 
ses  of  an  interconnect  with  an  insulation  envelope  have  also 
been  established.  Hence,  combined  electrochemical  and  heat- 
transfer  analyses  can  be  executed  through  an  iteration  process 
from  both  tools,  and  a  unified  temperature  distribution  with 
a  heat  loss  effect  can  be  obtained.  The  results  show  that  the 
convergent  rate  of  the  above  methodology  is  quite  efficient, 
and  that  the  temperature  patterns  are  compatible  with  coun¬ 
terparts  reported  in  the  literature. 

2.  Cell  performance  evaluation 

This  section  concerns  mainly  with  the  characteristics  of 
the  SOFC,  as  well  as  thermal,  flows  and  electric  parameters 
through  electrochemical  reactions.  Detailed  description  on 
the  computational  procedure  and  the  effects  of  various  pa¬ 
rameters  has  been  documented  elsewhere  [8]. 

2.7.  Electrochemical  computation 

Basically,  a  fuel  cell  is  an  electrochemical  reactor  that 
converts  chemical  energy  to  electrical  energy.  The  reactions 
at  the  anode  and  the  cathode  of  a  SOFC  that  operates  on 
hydrogen  and  oxygen  are  as  follows: 

at  anode  :  H2  +  O2-  —>  H2O  +  2e_  (1) 

at  cathode  :  ^C>2  +  2e_  — >  O2-  (2) 
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The  Nernst  potential  Fnst  from  chemical  reactions  can  be 
expressed  by: 


AG°  i  RT]n(Pa2P°?\ 
2  F  2 F  \  Ph2o  ) 


where  A G°  is  the  Gibbs  free  energy;  F  the  Faraday  constant; 
R  the  gas  constant;  and  Pu 2,  Pq2  and  Pn2o  are  the  partial 
pressures  of  H2,  O2  and  H2O,  respectively.  The  aforemen¬ 
tioned  potential  is  a  function  of  fuel  composition,  and  internal 
reforming  of  methane  will  be  considered  in  this  work. 

In  fact,  the  actual  operating  voltage,  Fop,  of  a  cell  will  be 
lower  than  the  theoretical  value  of  Fnst,  due  to  the  external 
connector  and  the  internal  resistance  of  the  PEN.  The  poten¬ 
tial  drop  within  the  cell  loop  considered  in  the  present  model 
includes:  (i)  ohmic  overpotential,  (ii)  concentration  overpo¬ 
tential  and  (iii)  activation  overpotential  [2,9-12].  Detailed 
formulations  of  these  overpotentials  have  been  documented 
previously  [8],  i.e., 
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Fig.  3.  Control  volume  of  SOFC  thermal  model. 


reactions,  i.e., 


Fop  —  Fast  Fpol  (4) 

Fpol  —  F0hm  T  Fact  T  FCOn  (5) 

The  chemical  reactions  between  the  hydrogen-rich  gas  and 
the  oxidant  are  the  major  mechanisms  that  govern  the  en¬ 
ergy  and  mass  conservations  in  the  cell.  For  a  natural  gas 
(NG)  fed  SOFC,  the  following  reactions  are  included  in  the 
electrochemical  model  and  are  also  shown  in  Fig.  1. 


ch4  +  h2o 

-+  3H2  +  CO 

(6) 

CO  +  h2o  - 

>  co2  +  h2 

(7) 

h2  +  io2  -> 

h2o 

(8) 

CO  + iQ2 - 

>  co2 

(9) 

Reaction  (6)  is  the  reforming  process  of  steam  with 
methane,  and  is  endothermic.  The  others  are  exothermic  re¬ 
actions.  The  heat  values  vary  with  temperature  and  can  be 
calculated  from  the  enthalpies  of  the  reactants  and  the  prod¬ 
ucts.  Reliable  and  comprehensive  thermodynamic  data  for 
this  work  were  obtained  from  the  JANAF  tables  [13].  These 
reactions  also  govern  the  mass  conservations  of  the  related 
chemical  species. 

2.2.  Thermal  model 

Based  on  energy  conservation  relations,  various  temper¬ 
atures  involved  in  the  present  model  can  be  calculated,  e.g., 
solid  temperature  Ts,  fuel  temperature  7f,  and  air  tempera¬ 
ture  Ta.  A  control  volume  is  shown  in  Fig.  3  and  includes 
fuel,  solid  and  air  units  for  the  thermal  model,  and  outlines 
the  relationships  among  heat  fluxes  and  source  terms.  Con¬ 
sidering  conduction  through  the  channel  wall  and  convection 
with  flow,  the  governing  equations  for  the  fluid  temperatures 
can  be  formulated  with  the  source  terms  from  the  chemical 


Y,  -  A{hf(Ts  -  7f)  =  Qf  (10) 

'  ox 

l 

Y  8(Q^P,Ta)dx  -  “  T a)  =  0a  (1 1) 

i 

In  Eqs.  (10)  and  (11),  C,  Cp,  A  and  h  represent  molar  flow 
rate,  specific  heat,  contact  area  and  heat-transfer  coefficient, 
respectively.  Moreover,  the  subscripts  /,  f,  s  and  a  indicate 
species  /,  fuel,  solid  and  air,  respectively. 

Similarly,  the  governing  equation  for  solid  temperature 
can  be  expressed  in  terms  of  conduction,  convection  and 
source  terms,  while  the  last  category  includes  chemical  en¬ 
ergy  from  the  reactions  and  the  Joule  heating  term  [10], 
namely: 

a  2rs  a2rs 

+  ky—fV, ev  -  Afhf(Ts  -  7f) 

-A&K{T,  -  TO  =  0S  (12) 

in  which  k  is  the  thermal  conductivity  of  the  solid,  Fcv  = 
AvAyAz,  and  the  source  term  is: 

Qs  —  Q\  (2a  —  f  (Fnst  F0hm  —  Fact  FCOn)  (13) 

To  circumvent  time-consuming  and  cost-intensive  com¬ 
putational  fluid  dynamics  efforts  to  model  convection  be¬ 
haviour,  empirical  correlations  for  heat-transfer  coefficients 
were  implemented  instead.  For  a  rectangular  channel,  the 
Nusselt  number,  Nu ,  changes  with  the  cross-sectional  aspect 
(width-to-height,  W:H)  ratio  [14].  In  preliminary  investiga¬ 
tions,  the  value  of  W\H  was  selected  to  be  around  2.5  and  3. 
Hence,  the  value  of  4  for  the  Nusselt  number  was  adopted 
to  calculate  the  heat-transfer  coefficients  for  the  fuel  and  air 
sides. 


130 


Y.-P.  Chyou  et  al.  /  Journal  of  Power  Sources  139  (2005)  126-140 


Table  1 

Inlet  and  outlet  conditions  of  SOFC 


Fuel 

Air 

Inlet 

Outlet 

Inlet 

Outlet 

Flow  rate  (mol  h  1  per  cell) 

1.20 

1.61 

20.00 

19.61 

Flow  rate  (mol  h_  1  per  stack) 

1.20  x  38 

1.61  x  38 

20.00  x  38 

19.61  x  38 

Composition  (%) 

ch4 

17.1 

0.11 

/ 

/ 

h2o 

49.34 

61.27 

/ 

/ 

h2 

26.26 

20.51 

/ 

/ 

CO 

2.94 

3.43 

/ 

/ 

co2 

4.36 

14.68 

/ 

/ 

n2 

/ 

/ 

79.00 

80.57 

o2 

/ 

/ 

21.00 

19.43 

S/C 

2.89 

/ 

o2/h2o 

/ 

7.09 

Temperature  (°C) 

Utilization(%) 

640 

66.53 

773 

640 

9.28 

763 

2.3.  Boundary  conditions  and  computation  procedures 

It  is  important  that  the  early  stage  of  research  and  devel¬ 
opment  is  compatible  with  the  state-of-the-art,  in  order  to 
understand  the  existing  status  of  technology  progress,  and  to 
benchmark  the  results  with  those  that  have  been  reported  in 
the  literature.  Hence,  this  work  has  adopted  the  ‘30%  pre¬ 
reforming  NG’  reported  by  Lehnert  et  al.  [15]  as  the  inlet 
conditions  of  the  SOFC,  for  which  the  compositions  are  listed 
in  Table  1 . 

To  simplify  the  analysis  model,  this  work  implemented 
an  idealized  and  bulk  approach  to  handle  the  boundary  con¬ 
ditions  in  the  reaction  zones.  It  has  been  reported  that  the 
pressure  loss  in  a  straight  type  of  flow  channel  is  usually  less 
than  5%;  hence,  this  effect  is  ignored  in  the  present  model. 
In  addition,  local  distribution  of  the  flow  field  is  ignored, 
and  mean  values  of  flow  speed  estimated  from  the  gas-flow 
rates  are  adopted  instead  in  the  reaction  zones.  Furthermore, 
the  adiabatic  temperature  distribution  is  first  calculated  from 
chemical  reactions,  and  then  the  effect  of  heat  loss  to  the 
environment  is  then  evaluated  (see  Section  3.3)  by  combin¬ 
ing  thermal  analysis  of  the  interconnect  via  a  finite  element 
model  (FEM). 

The  above  electrochemical  and  thermal  models  provide 
the  infrastructure  of  the  integrated  SOFC  analyses,  which  are 
implemented  into  an  indigenous  two-dimensional  (2D)  nu¬ 
merical  code.  The  solution  procedure  includes  several  steps, 
namely: 

(i)  specification  of  inlet  flow  conditions; 

(ii)  specification  of  geometry,  materials,  electric  parame¬ 
ters,  etc.; 

(iii)  estimation  of  initial  temperature  distributions  for  solid, 
fuel  and  air; 

(iv)  specification  of  current  density  to  calculate  the  operat¬ 
ing  voltage; 


(v)  utilization  of  the  Newton-Raphson  equation  to  per¬ 
form  electrochemical  computations; 

(vi)  evaluation  of  the  source  terms  and  heat-transfer  coef¬ 
ficients  for  the  energy  equations; 

(vii)  application  of  the  finite-volume  method  to  calculate 
temperature  distributions; 

(viii)  comparison  of  the  results  with  the  given  counterparts, 
and  implementation  of  iterative  processes  until  conver¬ 
gent  criteria  are  reached. 


3.  Thermal  analyses  of  the  interconnect 

The  interconnect  in  a  SOFC  stack  provides  several  ma¬ 
jor  functions;  namely:  (i)  flow  passages  for  working  me¬ 
dia  through  reaction  zones;  (ii)  an  electric  conductor  for  a 
serial  loop;  and  (iii)  structural  component  for  pile  assem¬ 
bly.  Due  to  the  above  functional  requirements,  it  is  desirable 
to  perform  thermo-mechanical  analyses  for  the  interconnect 
in  order  to  ensure  its  durability  and  mechanical  integrity. 
The  early  stage  of  planning  for  this  issue  has  been  focused 
on  establishing  a  finite  element  model  and  a  computational 
procedure  [16]. 

3.1.  Finite  element  model 

As  shown  in  Fig.  2,  the  present  interconnect  design  is  a 
rectangular  plate  with  a  length  of  136  mm,  a  width  of  95  mm 
and  a  thickness  of  3  mm.  Flow  channels  are  embedded  in  the 
central  region  of  the  plate.  These  provide  passages  for  fuel 
and  oxidant  and  are  subjected  to  a  thermal  load.  To  reduce 
heat  loss  to  the  environment  and  the  outer  temperature  of  the 
SOFC,  the  unit  is  enclosed  by  two-layered  insulation  ma¬ 
terials.  The  basic  configuration  consisted  of  a  10-mm  inner 
layer  of  ceramic  fibre  and  a  7 5 -mm  ceramic  plate  as  the  outer 
insulation. 
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Fig.  4.  Finite  element  model  of  interconnect  with  thermal  insulations. 

Based  on  the  above  dimensions  of  the  interconnect  and  the 
insulation  layers,  a  finite  element  model  has  been  established 
by  using  the  generic  computer-aided  engineering  (CAE)  code 
of  ANSYS  [17].  This  solid  model  has  been  built  up  with 
the  SOLID70  element,  and  includes  9760  and  574  elements 
for  the  interconnect  and  the  insulation  layers,  respectively. 
The  finite  element  model  of  the  interconnect  with  thermal 
insulation  layers  is  presented  in  Fig.  4,  and  the  3D  mesh  of 
the  interconnect  is  shown  in  Fig.  5. 

The  provisional  interconnect  was  made  of  Incoloy-800. 
The  material  properties  are  documented  in  [16]  and  the  ther¬ 
mal  conductivity  at  800  °C  (i.e.,  k\  =  25  Wm-1  °C_1)  was 


Fig.  5.  3D  finite  element  model  of  interconnect. 


adopted  in  the  present  analysis.  In  addition,  constant  val¬ 
ues  of  thermal  conductivity  were  used  for  the  insulation,  i.e., 
k,2  =  0.198Wm_1  °C_1  for  the  inner  ceramic  fibre  and  k 3 
=  0.215  Wm-1  °C-1  for  the  outer  ceramic  plate.  The  sur¬ 
rounding  temperature  was  assumed  to  be  40  °C. 

3.2.  Two-dimensional  heat-transfer  analysis 

Based  on  the  calculated  results  in  Section  2,  the  tem¬ 
perature  distribution  has  a  2D  nature  in  the  reaction  zones. 
Hence,  the  corresponding  heat-transfer  behaviour  in  the 
interconnect  will  exhibit  2D  features  as  well.  In  practice,  the 
stack  is  enveloped  by  insulation  materials  of  low  thermal 
conductivity  to  diminish  the  heat  loss  to  the  environment  and 
the  temperature  at  the  outer  edges.  Accordingly,  the  thermal- 
transport  mechanism  within  the  solid  parts  is  dominated  by 
heat  conduction  and  this  substantially  relieves  the  demand 
of  CFD  efforts.  At  the  outer  boundary  of  the  insulation 
envelop,  heat  transfer  is  governed  by  natural  convection  to 
the  surroundings.  Due  to  relatively  lower  temperatures  in 
this  region,  the  effect  of  local  gradients  is  minimal;  hence, 
heat-transfer  coefficients  are  estimated  by  taking  the  average 
temperature.  This  simplified  approach  should  pose  only 
a  minor  influence  on  the  heat-transfer  behaviour  in  the 
cell  core. 

In  summary,  heat-transfer  phenomena  at  the  interconnect 
can  be  simplified  to  ‘quasi-conduction’  analyses,  for  which 
solutions  are  directly  calculated  by  utilizing  thermal  elements 
embedded  in  the  CAE  code.  The  benefits  of  this  approach 
are  two-fold.  On  one  hand,  the  temperature  distribution 
at  the  interconnect  can  be  determined  and  the  necessary 
boundary  conditions  for  subsequent  thermal  stress  analysis 
can  be  obtained.  On  the  other  hand,  the  amount  of  heat 
loss  can  be  evaluated  and  fed  back  to  the  ECM  for  further 
evaluation  and  the  effect  of  insulation  thickness  can  also  be 
studied. 

3.3.  Iteration  between  FEM  and  ECM 

The  integrated  thermal  analyses  of  a  planar  SOFC  were 
divided  into  two  stages.  First,  the  temperature  distribution 
at  the  reaction  zone  of  the  interconnect  was  evaluated  via 
ECM  under  adiabatic  conditions,  and  was  utilized  as  the  ther¬ 
mal  boundary  for  the  FEM.  Second,  2D  heat- transfer  analy¬ 
sis  was  performed  to  obtain  the  overall  thermal  parameters. 
These  included  a  heat-loss  term  that  was  later  fed  back  to  the 
cell  performance  evaluation  for  a  more  realistic  considera¬ 
tion.  Through  several  iterative  calculation  processes,  a  uni¬ 
fied  temperature  distribution  could  be  obtained.  Details  of  the 
calculation  procedure  and  related  parametric  data  have  been 
reported  elsewhere  [6]. 

The  flow  chart  of  the  solution  procedure  is  presented 
in  Fig.  6.  At  the  first  step,  the  adiabatic  boundary  condi¬ 
tion  is  implemented  for  the  electrochemical  calculations  in 
the  reaction  zone  until  a  convergent  solution  is  achieved. 
Then,  the  FEA  is  activated.  The  primary  adiabatic  solu- 
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Fig.  6.  Flowchart  of  SOFC  electrochemical  model. 


SOFC  will  be  outlined,  and  a  model  for  design  analysis  will 
also  be  illustrated. 


4.1.  Bulk  heat-transfer  model 


The  approach  here  is  to  implement  the  logics  of  theo¬ 
retical  analyses,  i.e.,  first  to  clarify  the  relationship  among 
various  heat-transfer  mechanisms  via  related  empirical  cor¬ 
relations,  then  to  perform  one-dimensional  (ID)  analytical 
solutions  with  bulk  values.  In  the  overall  model,  the  heat 
source  originates  from  the  hot  gas  temperature  in  the  reac¬ 
tion  zones;  while  the  heat  flux  is  conveyed  to  the  channel 
walls  through  forced  convection,  and  to  the  edges  of  the  inter¬ 
connect  via  conduction.  Heat  transfer  through  the  insulation 
layers  around  the  stack  is  also  governed  by  conduction,  and 
to  the  environment  by  natural  convection. 

Hence,  the  governing  equation  of  heat  flux  can  be  ex¬ 
pressed  as  follows: 


Twi)  —  k\ 


Tw2 

AXi 


T\v2  Tw  3 
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ax3 


ho(Tw4  Too)- 


(14) 


where  the  subscripts  1 ,  2  and  3  for  k  and  AX  represent  the 
interconnect,  the  inner  insulation  and  the  outer  layer,  respec¬ 
tively;  Ti  and  are  the  fluid  temperature  inside  the  channels 
and  the  outer  surroundings,  respectively;  and  Tw\  to  Tw 4  are 
the  wall  temperatures  at  the  edges  of  solid  parts. 

Introducing  the  concept  of  total  heat  resistance,  Eq.  (14) 
can  be  rewritten  as  follows: 
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tion  is  utilized  as  the  boundary  condition  for  heat-transfer 
computation  in  the  second  step.  Iterations  between  electro¬ 
chemical  and  heat-transfer  analyses  are  then  performed  until 
convergence. 


where  U  is  the  total  heat-transfer  coefficient;  A  the  total  heat- 
transfer  area;  R  the  total  heat  resistance;  and  AT  =  Ti  — 
Tq O’ 

4.2.  Correlations  of  forced  convection 


4.  Model  for  design  analysis 

The  models  described  above  provide  the  basis  for  in¬ 
tegrated  thermal  analyses  of  a  planar  SOFC.  It  would  be 
favourable,  however,  to  develop  simplified  design  analysis 
tools  for  component  prototyping,  prior  to  employing  compre¬ 
hensive  finite  element  analyses.  For  example,  the  intercon¬ 
nect  is  an  important  structural  component  of  SOFC  stacks  and 
its  thermal  characteristics  are  essential  for  subsequent  analy¬ 
ses  of  mechanical  integrity  in  long-term  operation.  Based  on 
the  geometry  and  sizing  of  the  flow  channels  on  the  intercon¬ 
nect,  some  crucial  features  of  heat-transfer  behaviour  in  the 


The  reaction  zones  of  a  planar  SOFC  consist  of  a  series 
of  parallel  flow  channels.  Thus,  the  heat-transfer  behaviour 
from  the  internal  source  to  the  environment  may  be  regarded 
as  laminar  flow  over  a  flat  plate,  from  a  global  viewpoint. 
Forced  convection  phenomena  can  be  simulated  with  empir¬ 
ical  correlations  documented  in  the  literature  [18]. 

For  external  flow  over  a  flat  plate,  the  mean  Nusselt  num¬ 
ber  that  governs  the  heat-transfer  behaviour  normal  to  the 
plate  can  be  correlated  to  the  Reynolds  number,  Re ,  and  the 
Prandtl  number,  Pr ,  in  the  following  expression: 

Num  i  =  —  =  0.664 Rel/1Prl/3  (16) 

ki 


Y.-P.  Chyou  et  al.  /  Journal  of  Power  Sources  139  (2005)  126-140 


133 


where 

uL 

Re  =  — 
v 


(17) 


In  the  above  equations,  ki,  Pr  and  v  denote  the  ther¬ 
mal  conductivity,  Prandtl  number  and  kinematic  viscosity 
of  the  fluid,  respectively.  In  addition,  L  is  the  characteristic 
length  of  the  flow  field  and  u  is  the  averaged  velocity.  Hence, 
the  heat-transfer  coefficient  for  forced  convection  hi  can  be 
calculated. 

The  Nusselt  number  will  approach  a  constant  value  for 
internal  flow,  provided  that  a  fully-developed  flow  condition 
is  achieved.  Two  types  of  flow  configuration  will  be  illus¬ 
trated  in  the  following.  The  first  case  (Nu\)  is  concerned  with 
flow  confined  by  two  plates,  one  of  which  is  adiabatic.  The 
second  example  (Nu2)  deals  with  heat  transfer  through  the 
channel  walls,  i.e.,  confined  flow.  The  former  governs  the 
heat- transfer  behaviour  normal  to  the  plate,  i.e.,  similar  to  that 
given  by  Eq.  (16),  which  should  be  appropriate  in  discussing 
the  situation  at  the  end  plates;  while  the  latter  deals  with  ‘all- 
around’  heat  flux  through  a  confined  channel,  which  might 
be  utilized  to  simulate  the  lateral  heat  flux  at  the  periphery  of 
the  cell.  As  mentioned  earlier,  Nu2  for  a  rectangular  channel 
changes  with  the  cross-sectional  W:H  ratio;  nevertheless,  the 
values  listed  in  the  literature  deviate  only  slightly.  Hence,  the 
value  of  Nu2  given  in  the  textbook  written  by  Ozisik  [18]  will 
be  adopted  here,  for  easier  reference  purposes,  for  the  design 
analysis  model. 

The  correlations  are  listed  here  for  further  discussion. 
Since  the  flow  rate  of  air  is  substantially  larger  than  that  of 
the  fuel,  heat-transfer  coefficients  in  the  air  channels  will  be 
implemented  in  the  ID  simplified  analyses,  as  a  more  con¬ 
servative  approach. 


Nu  i  =  5.385  for  flow  between  plates  (one  side  adiabatic) 


(18) 


Nu2  =  3.68  for  flow  confined  in  channel  (19) 

4.3.  Correlations  of  natural  convection 

Natural  convection  phenomenon  can  be  simulated  with 
the  following  correlation. 

Num,0  =  =  0.5 14 (GrL  Pr)0'25  (20) 
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Fig.  7.  Comparison  of  cell  performance. 


5.  Results  and  discussion 

5.1.  Cell  performance 

5.1.1.  Influence  of  flow  arrangement 

There  are  usually  three  types  of  flow  arrangement  for  a 
planar  SOFC,  i.e.,  co-flow,  counter-flow  and  cross-flow.  Cell 
characteristics  will  be  affected  by  different  flow  patterns.  A 
comparison  of  power  density  and  temperature  for  various 
flow  arrangements  under  the  same  inlet  conditions  is  given 
in  Fig.  7.  It  is  seen  that  the  counter-flow  arrangement  results 
in  the  highest  power  density  and  temperature,  which  is  the 
primary  configuration  studied  in  this  paper. 

Further  details  of  temperature  patterns  in  the  active  area  of 
a  cell  with  various  flow  arrangements  are  revealed  in  Fig.  8. 
For  simplicity,  an  adiabatic  condition  is  applied  at  the  bound¬ 
aries  of  the  computation  domain,  which  might  affect  the  tem¬ 
perature  contours  near  the  periphery  but  would  not  distort  the 
generic  feature  in  the  core  area.  It  is  seen  that  the  cross-flow 
configuration  generated  a  hot  island  near  the  fuel  inlet  and  the 
air  outlet  in  the  reaction  zone,  as  shown  in  Fig.  8a.  More  uni¬ 
form  temperature  distributions  can  be  obtained  with  parallel- 
flow  configurations.  The  temperature  pattern  of  the  counter¬ 
flow  configuration  is  shown  in  Fig.  8b,  while  the  counterpart 
of  the  co-flow  case  is  given  in  Fig.  8c.  These  two  cases  dis¬ 
play  a  similar  pattern,  i.e.,  the  temperature  increases  along  the 
airflow  direction;  but  differ  in  terms  of  temperature  gradient 
and  level.  In  summary,  the  outcome  obtained  from  the  ECM 
in  this  work  is  qualitatively  consistent  with  that  reported  in 
the  literature,  e.g.  [5]. 
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In  the  above  equation,  k0  is  the  thermal  conductivity  of  the 
ambient  atmosphere;  hence,  the  heat- transfer  coefficient  for 
natural  convection  hQ  can  be  calculated. 


5.1.2.  Electrochemical  characteristics  of  SOFC 

One  of  the  major  considerations  in  this  work  has  been 
the  utilization  of  metal  rather  than  ceramics  as  the  material 
for  the  interconnect,  in  order  to  reduce  cost  as  well  as  to 
enhance  mechanical  integrity  and  durability.  Hence,  it  has 
been  set  in  the  project  specifications  that  the  maximum  metal 
temperature  should  be  kept  below  or  around  800  °C.  Through 
several  trial  calculations  with  the  ECM  described  in  Section 
2,  the  inlet  temperature  of  640  °C  roughly  complies  with  the 
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Fig.  8.  (a)  Temperature  pattern  of  cross-flow  cell  with  adiabatic  boundaries,  (b)  Temperature  pattern  of  counter-flow  cell  with  adiabatic  boundaries,  (c) 
Temperature  pattern  of  co-flow  cell  with  adiabatic  boundaries. 


project  specifications.  The  resulting  steady-state  average  and 
maximum  solid  temperatures  from  the  model  computation 
are  773  and  802  °C,  respectively.  Fuel  utilization  is  66.5%, 
while  the  air  consumption  rate  is  at  a  minimal  level  of  below 
10%  due  to  a  relatively  high  value  of  the  oxygen  to  water 
ratio  (X  ~  7).  The  input  and  output  data  are  summarized  in 
Table  1  for  further  reference. 

The  corresponding  temperature  distribution  in  the  reaction 
zone  is  presented  in  Fig.  9.  The  distributions  of  fluid  temper¬ 
atures  exhibit  similar  patterns,  while  the  levels  are  slightly 
different.  Near  the  fuel  inlet,  the  maximum  values  of  Ts,  7f 
and  Ta  are  802.4,  806.3  and  800.6  °C,  respectively.  The  value 
of  ra  at  the  air  inlet  is,  however,  somewhat  lower  than  the  oth¬ 
ers;  the  minimum  values  of  Ts,  7f  and  Ta  are  7 16.9, 721 .6  and 
680.3  °C,  respectively.  It  should  be  noted  that  the  tempera¬ 
ture  pattern  shown  in  Fig.  9  is  the  outcome  of  the  integrated 


electrochemical  and  heat  transfer  analyses  illustrated  in  this 
work,  which  relieve  the  constraint  of  an  adiabatic  boundary 
condition  at  the  edges  of  the  reaction  zone.  On  comparing  the 
results  in  Fig.  9  with  those  in  Fig.  8b,  it  is  observed  that  the 
adiabatic  boundaries  impose  a  certain  influence  on  the  tem¬ 
perature  distribution  at  the  edges  of  the  active  area.  By  con¬ 
trast,  however,  the  central  region  remains  almost  unaffected, 
as  far  as  the  contour  patterns  are  concerned.  This  outcome 
confirms  the  statement  mentioned  in  the  previous  section.  In 
addition,  it  will  be  shown  later  in  this  paper  that  the  temper¬ 
ature  level  will  be  somewhat  lower  in  the  former  case,  due  to 
the  overall  heat  loss. 

The  overall  cell  characteristics  are  presented  in  Fig.  10. 
As  reported  in  various  studies,  the  value  of  0.7  V  has 
been  selected  as  the  operating  voltage  of  a  single  cell  for 
comparison  purposes.  The  data  yield  630  mA  cm-2  and 
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Fig.  9.  Temperature  distribution  in  cell. 


440  mW  cm  2  for  the  current  density  and  the  power  density, 
respectively. 

5.2.  Thermal  characteristics  of  the  interconnect 

Due  to  the  symmetric  nature  of  the  present  design,  only 
half  of  the  interconnect  has  been  simulated  in  the  analysis 
model;  the  other  half  can  be  presented  as  the  mirror  image. 
As  mentioned  above,  a  unified  temperature  distribution  from 
both  the  electrochemical  and  finite  element  analyses  can  be 
obtained  after  several  iterations,  which  will  be  shown  in  this 
section.  Global  temperature  data  of  the  finite  element  model 
for  the  interconnect  with  a  two-layered  insulation  envelope 
are  given  in  Fig.  11.  It  is  seen  that  steep  temperature  gra¬ 
dients  exist  across  the  ceramic  layers,  which  confirms  the 
function  of  the  insulation.  In  general,  the  outer  temperature 
of  the  insulation  envelope  adjacent  to  the  atmosphere  can  be 
reduced  to  values  below  200  °C,  depending  on  the  material 
properties  and  thickness  of  the  multi-layers.  Compared  with 
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Fig.  10.  SOFC  characteristics. 


the  SOFC  operating  temperature  of  around  800  °C  in  the  ac¬ 
tive  area,  lower  temperatures  of  this  level  can  substantially 
mitigate  the  amount  of  heat  loss  to  the  environment.  Never¬ 
theless,  layout  or  optimization  of  stack  insulation  is  subject 
to  engineering  practice. 

Further  insights  of  the  thermal  characteristics  can  be  re¬ 
vealed  by  examining  the  FEA  results  for  the  interconnect 
itself,  as  presented  in  Fig.  12.  In  addition,  the  heat  flux  distri¬ 
bution  at  the  periphery  of  the  active  area  as  well  as  the  value 
of  the  heat  loss  per  unit  cell  can  be  calculated  by  the  FEA. 
Information  about  the  heat  flux  distribution  of  the  compu¬ 
tation  domain,  i.e.,  upper  half  of  the  interconnect  shown  in 
Fig.  12,  is  presented  in  Table  2.  It  is  observed  that  heat  flux 
normal  to  the  flow  channels,  qy ,  varies  from  positive  values 
near  the  fuel  inlet  to  negative  values  towards  the  air  inlet. 
This  phenomenon  shows  a  generic  pattern  in  the  active  area 
with  parallel-flow  arrangements,  i.e.,  the  lateral  heat  flux  is 
outwards  near  the  fuel  inlet  side,  but  turns  inwards  on  ap¬ 
proaching  the  air  inlet  region.  The  heat  flux  parallel  to  the 
flow  channels,  qx ,  is  negative  at  the  fuel  inlet  but  positive  at 
the  opposite  edge,  which  means  that  both  the  fuel  and  the 
air  flow  encounter  upward  temperature  gradients  along  their 
individual  flow  direction  for  the  counter-flow  configuration. 
In  summary,  the  results  indicate  two  major  features:  (i)  heat 
generation  due  to  chemical  reactions  is  most  intense  near  the 
fuel  inlet;  and  (ii)  airflow  is  most  effective  at  cooling  near  the 
air  inlet. 

5.3.  Computational  aspects 

As  described  in  Section  3.1,  two-layered  insulation  mate¬ 
rials,  with  an  aim  to  reduce  the  heat  loss  to  the  environment, 
enclosed  the  SOFC  unit.  Three  different  insulation  envelopes 
have  been  implemented  in  this  work  to  investigate  the  heat 
transfer  behaviour  of  a  metallic  interconnect  with  a  ceramic 
enclosure.  The  inner  layer  of  ceramic  fibre  was  kept  the  same 
thickness  of  10  mm,  while  three  different  values  of  50,  75 
and  100  mm  were  selected  for  the  outer  plate.  These  cases 
are  identified  here  as  A,  B  and  C,  respectively.  Furthermore, 
it  should  be  noted  that  case  B  is  the  basic  configuration  illus¬ 
trated  in  Fig.  12.  The  temperature  differences  (AT)  between 
the  results  from  ECM  and  FEA  for  various  insulation  ar¬ 
rangements  are  presented  in  Fig.  13a-c.  It  is  seen  that  the 
value  of  AT  falls  to  less  than  1  °C  from  a  bulk  value  around 
1000  °K  on  the  interconnect,  after  three  iterations.  This  out¬ 
come  indicates  that  an  accuracy  of  better  than  0.1%  has  been 
achieved  for  all  three  cases.  Hence,  it  has  been  confirmed  that 
the  convergent  rate  of  the  iteration  process  is  very  efficient. 

For  reference  purposes,  the  computational  resources  re¬ 
quired  for  the  analyses  are  given  here.  The  electrochemical 
calculations  were  performed  on  a  Sun  Ultra2  workstation, 
with  167  MHz  CPU  and  128  MB  RAM,  and  the  CPU  time 
was  around  1.5  h  for  each  primary  analysis.  For  subsequent 
iterations,  much  less  computational  efforts  were  required.  For 
example,  the  CPU  time  for  the  first  three  iterations  turned  out 
to  be  around  3,  2.5  and  2  min,  respectively.  Alternatively,  the 
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STEP=1 
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PowerGraphics 
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Fig.  11.  Temperature  distribution  of  interconnect  including  insulation. 


Fig.  12.  Temperature  distribution  on  interconnect. 


Table  2 


Heat  flux  on  interconnect  boundary 


Item 

1  2 

3 

4 

5 

6 

7 

8 

9 

10 

11 

12 

13 

14 

At  upper  periphery  (kW  m-2) 
x(mm)  1.93  7.85 

qy  3.37  5.13 

13.35 

6.5 

18.85 

7.54 

24.35 

8.1 

29.85 

8.05 

35.35 

7.28 

40.85 

5.74 

46.35 

3.51 

51.85 

0.8 

57.35 

-1.96 

62.85 

-4.16 

68.35 

-4.89 

74.28 

-2.64 

On  the  fuel  inlet  (left)  side  (kW 

m-2) 

On  the  air  inlet  (right)  side  (kW  m  2) 

8a  9a 

10a 

lla 

12a 

13a 

14a 

8a 

9a 

10a 

lla 

12a 

13a 

14a 

y  (mm) 

4x 

40.85  46.35 

-2.91  -4.33 

51.85 

-4.47 

57.35 

-4.35 

62.85 

-4.15 

68.35 

-3.79 

74.28 

-3.46 

40.85 

3.66 

46.35 

0.97 

51.85 

0.54 

57.35 

0.56 

62.85 

0.57 

68.35 

0.55 

74.28 

0.37 

a  Item. 
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1  st  layer:  I  Omm,  2nd  layer:50mm 


1st  layer:  10mm,  2nd  layer:  75mm 


(b)  Iteration  number 

—♦—min.  -»-max. 


1  st  layer:  I  Omm,  2nd  layer:  1 00mm 


(c)  Iteration  number 

min.— ■—  max. 

Fig.  13.  (a)  Convergent  rate  of  iterations:  case  A.  (b)  Convergent  rate  of 
iterations:  case  B. 

platform  for  FEA  was  a  PENTIUM  IV,  32-bit  CPU  3.06  GHz, 
RAM  768  MB,  and  each  computation  needed  a  CPU  time  of 
slightly  more  than  2  min. 

5.4.  Parameter  sensitivities 

As  described  in  Section  4,  simplified  design  analysis  tools 
would  be  favourable  in  the  prototyping  process,  i.e.,  for  pro¬ 
viding  a  preliminary  database  of  related  components  as  well 
as  for  efficient  parametric  studies.  Hence,  the  sensitivities  of 
various  thermal  and  flow  parameters  would  be  of  major  con¬ 
cern  in  identifying  the  predominant  governing  mechanisms. 

5.4.1.  Heat  resistance 

From  Eq.  (15),  it  is  noted  that  the  heat  transfer  behaviour 
from  cell  to  environment  is  governed  mainly  by  the  following 
eight  parameters:  hi ,  k\,  &2,  &3,  h0 ,  Ax\,  Ax2  and  AV3.  The 
total  heat  resistance,  ^  R ,  consists  of  five  entities:  (1)  R\  = 
Ax\/k\ ,  (2)  R2  =  Ax2 /&2>  (3)  R3  =  AV3/&3,  (4)  Ri  =  l/hi  and 


(5)  R0  =  1  /h0.  The  first  three  items  represent  the  resistances 
due  to  conduction  in  various  materials,  while  the  latter  two 
are  the  counterparts  from  convection.  The  basic  configura¬ 
tion  mentioned  above,  i.e.,  case  B,  has  been  adopted  here 
as  a  paradigm  to  illustrate  the  relative  influence  of  the  var¬ 
ious  heat  resistances.  For  appropriate  boundary  conditions, 
the  temperature  of  hot  gas  within  the  cell  and  surrounding  at¬ 
mosphere  is  assumed  to  be  800  and  40  °C,  respectively.  In  the 
reference  case  for  lateral  heat  flux  at  the  periphery  of  the  cell, 
Eq.  (19)  was  implemented  to  evaluate  the  heat-transfer  coef¬ 
ficient  with  internal  forced  convection,  i.e.,  hi.  Nevertheless, 
other  correlations  will  be  discussed  as  well  in  a  later  section 
of  this  paper.  In  addition,  the  thickness  between  the  outer¬ 
most  flow  passage  to  the  edge  of  interconnect  was  selected 
as  the  value  for  the  channel  wall,  as  a  conservative  measure 
for  the  resistance  to  heat  conduction  in  the  metal. 

All  the  aforementioned  geometric,  material,  thermal  and 
flow  parameters  are  summarized  in  Table  3  for  easy  reference. 
Calculated  values  of  the  heat  resistances  are  also  given.  It  is 
seen  that  the  heat  resistance  of  the  outer  insulation  layer, 
R3 ,  assumes  the  highest  value  for  the  present  configuration. 
Hence,  the  predominant  parameters  to  affect  the  heat-transfer 
behaviour  should  be  the  thickness  and  material  property  of  the 
ceramic  plate  used  for  insulation.  This  outcome  also  confirms 
the  adequacy  of  applying  an  insulation  envelope  to  enclose 
the  SOFC  unit,  for  reducing  heat  loss  to  the  environment. 

5.4.2.  ID  heat-transfer  results 

Parametric  studies  with  the  ID  bulk  model  for  design  anal¬ 
ysis  described  in  Section  4  have  been  performed  to  investigate 
the  influence  of  insulation  layer  thickness  as  well  as  heat- 
transfer  coefficient.  The  test  matrix  of  ID  analyses  covered 
the  combination  of  four  different  values  of  A3C3  and  various 
correlations  of  hi.  The  data  obtained  from  the  simplified  ap¬ 
proach  were  then  compared  with  the  counterparts  of  FEM 
reported  in  Section  5.2,  and  empirical  correlations  could  be 
deduced  as  the  ‘accommodation  (rule-of-thumb)  factors’  for 
rapid  prototyping  purposes. 

Referring  to  Fig.  11,  the  global  2D  temperature  distribu¬ 
tion  shows  a  pattern  similar  to  that  of  a  ‘rectangular  duct’ 
enclosed  by  insulation  materials.  Hence,  the  concept  of  an 


Table  3 

Parameters  of  one-dimensional  heat  transfer 


Geometry  (mm) 

Axi 

AX2  AX3 

L  H  W 

9 

10  75 

76  1  3 

Thermal  conductivity  (W  m  1  °C  l) 

Heat  transfer  coefficients 

(Wm-1  “C”1) 

h  k  2 

£3 

hi  hQ 

25  0.198 

0.215 

154.6  9.88 

Heat  transfer  resistance  (m2  °C  W  1 ) 

R\  =  Ax\/k\ 

R2  =  AX2/&2  R3  = 

Ax3/&3  Rf  =  1/hi  Ro  =  VhQ 

0.00017 

0.051  0.349  0.0065  0.101 
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Table  4 

Results  of  one-dimensional  heat  transfer 


AX3  (cm) 

Tv  1  (°C) 

Tv 2  (°C) 

Tv 3  (° 

C)  Tv4  (°C) 

q  (kW  m  2) 

hi 

=  11.85  W  m-2  °C 

1  (plate),  Num  = 

11.5,  L  =  76  mm 

2.5 

618.2 

617.4 

508.6 

258.1 

2.15 

5.0 

663.3 

662.7 

580.8 

204.0 

1.62 

7.5 

690.4 

690.0 

624.4 

171.4 

1.30 

10 

708.6 

708.2 

653.5 

149.6 

1.08 

hi 

=  4.46  Wm-2  °C“ 

1  (adiabatic 

:,  plate),  Nu\  =  5.385,  L  =  76  mm 

2.5 

454.1 

453.5 

375.6 

196.2 

1.54 

5.0 

520.1 

519.7 

456.6 

166.3 

1.25 

7.5 

565.0 

564.6 

511.7 

146.1 

1.05 

10 

597.5 

597.2 

551.5 

131.4 

0.90 

hi 

=  154.6Wm-2  °C 

1  (channel),  Nu2 

=  3.68,  Dh  = 

1.5  mm 

2.5 

782.1 

781.1 

641.5 

319.9 

2.77 

5.0 

787.4 

786.7 

688.6 

236.7 

1.94 

7.5 

790.3 

789.8 

714.1 

191.6 

1.50 

10 

792.1 

791.7 

730.1 

163.3 

1.22 

Thermal 

boundary 

conditions: 

Tmax 

=  800 °C, 

T0  =  40  °C,  h0  = 

9.88  W m-2  °C_1. 


‘outer  critical  insulation  radius  (roc)’  might  apply  for  check¬ 
ing  the  suitability  of  the  insulation  thickness,  provided  that  an 
effective  diameter  of  the  ‘rectangular  duct’  can  be  deduced 
and  can  replace  the  diameter  of  a  circular  tube  in  a  conven¬ 
tional  approach.  As  given  in  textbooks,  e.g.  [18],  the  roc  for 
a  tube  is  a  function  of  the  thermal  conductivity  of  the  insula¬ 
tion  material,  ki ,  and  the  heat-transfer  coefficient  at  the  outer 
boundary,  hQ ,  and  can  be  expressed  as  follows: 


1D-Plate 

ID-Channel 


1  D-Adiabatic 


— 2D  FEM 


Ax3  [mm] 


Fig.  14.  (a)  Heat  flux  vs.  insulation  thickness,  (b)  Relationship  of  lateral 
heat  flux. 


r  oc 


ki_ 

h0 


(23) 


On  inserting  the  values  of  k 3  and  h0  into  Eq.  (23),  roc  is  found 
to  be  around  22  mm.  This  outcome  indicates  that  the  critical 
insulation  diameter  (2  x  roc)  is  smaller  than  the  dimensions  of 
the  interconnect.  Therefore,  the  physical  configuration  in  this 
work  has  already  gone  beyond  the  regime  of  critical  insula¬ 
tion,  and  heat  loss  should  diminish  with  increasing  thickness 
of  the  insulation. 

The  results  obtained  from  the  ID  bulk  heat- transfer  anal¬ 
yses  are  summarized  in  Table  4,  in  which  wall  temperatures 
and  heat  flux  are  listed  against  AV3  and  hi.  For  convenience, 
the  related  data  extracted  from  cases  A,  B  and  C  mentioned 
above  are  also  presented  in  Table  5  for  comparison  with  the 
counterparts  from  ID  analysis  for  lateral  heat  flux.  Some  se¬ 
lected  data  from  Tables  4  and  5  are  presented  in  graphical 


Table  5 


Summary  of  2D  FEA  results 


Ax3  (cm) 

Tv 4  (°C) 

q  (kW  m  2), 
top 

q  (kW  m  2), 
right 

q  (kW  m  2), 
left 

Two-dimensional  heat  transfer 

5.0  177.5  3.36 

1.67 

-4.70 

7.5 

144.7 

2.89 

0.97 

-4.04 

10 

121.7 

2.62 

0.46 

-3.63 

ECM:  hi  =  203.5  Wm-2  °C_1,  kf  =  0.0763  Wm-1  °C“1,  Nu  =  4,  Dh  = 
1.5 mm.  FEA:  kx  =  25  Wm-1  0CT1,  Ajc2  =  1  cm,  k2  =  0.198  Wm"1  °C"1, 
k3  =0.215  Wm”1  °C_  1 ,  hQ  =  9.88  W m— 2  °C— 1 ,  T0  =  40°C. 


forms  for  easier  illustration  and  interpretation.  In  engineer¬ 
ing  practice,  the  heat  loss  to  the  environment  and  the  outer 
material  temperature  of  the  insulation  envelope  are  of  major 
concern  and  are  shown  in  Figs.  14  and  15,  respectively. 


(a) 


ID-Plate 
ID-Channel 


Q 

<M 


I 

0.9 

0.8 

0.7 
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50 


(b) 


Ax3  [  mm] 


1  D_Adiabatic 
2D  FEM 


-A —  1  D-Channel 


A 


75 

Ax 3  [mm] 


100 


Fig.  15.  (a)  Outer  material  temperature  vs.  insulation  thickness,  (b)  Rela¬ 
tionship  of  material  temperature. 
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The  heat  flux  values  of  ID  calculations  based  on  three  dif¬ 
ferent  correlations  as  well  as  the  lateral  counterpart  of  the  2D 
FEA  are  shown  in  Fig.  14a.  In  general,  the  heat  flux  decreases 
with  increasing  thickness  of  the  insulation  layer,  as  expected. 
In  addition,  the  variation  of  the  heat-transfer  coefficient  for 
internal  forced  convection  is  somewhat  neutralized  by  the 
larger  heat  resistances  of  the  insulation  materials;  hence,  the 
values  of  heat  flux  exhibit  relatively  minor  changes  among 
various  ID  models.  It  should  be  noted  that  heat  transfer  over  a 
flat  plate,  either  in  an  open  field  or  confined,  mainly  concerns 
the  heat  flux  normal  to  the  plate,  which  would  be  suitable  for 
simulating  the  thermal  behaviour  on  the  end  plates.  The  rea¬ 
son  for  presenting  the  results  based  on  the  correlations  of 
Numj  and  Nu\  in  Fig.  14a  is  two-fold:  on  the  one  hand,  it 
shows  the  generic  feature  and  capability  of  the  model  illus¬ 
trated  in  Section  4;  on  the  other  hand,  these  data  are  required 
in  later  work  to  investigate  the  heat-transfer  behaviour  of 
the  whole  stack  rather  than  a  single  cell.  Geometrically,  the 
flow  channels  on  the  interconnect  are  rectangular  ducts.  At 
present,  this  paper  focuses  on  the  heat-transfer  phenomenon 
at  the  periphery  of  an  SOFC;  hence,  the  correlation  of  chan¬ 
nel  flow  (i.e.,  Nu2)  would  be  more  appropriate,  although  not 
identical  to  the  present  physical  configuration,  in  establishing 
the  simplified  model  for  design  analysis. 

As  illustrated  in  Fig.  12,  the  thermal  characteristics  of 
planar  cells  exhibit  a  2D  nature.  Therefore,  some  ‘accom¬ 
modation  factors’  would  be  feasible  to  correlate  the  present 
ID  results  with  the  counterparts  of  2D  FEA.  The  relation¬ 
ship  between  the  lateral  heat  flux  and  the  thickness  of  the 
insulation  material  is  presented  in  Fig.  14b.  The  accommo¬ 
dation  factor  of  heat  flux,  AF^  =  g(lD)/g(2D),  is  found  to  vary 
quasi-linearly  within  a  range  between  0.58  and  0.48,  forAv3 
increasing  from  50  to  100  mm  in  the  present  investigation. 
This  feature  is  quite  promising  from  design  analysis  view¬ 
points.  More  parametric  studies  could  be  performed  in  the 
future  to  extend  the  comprehensiveness  of  the  correlations, 
which  would  be  beneficial  to  rapid  prototyping  processes. 

Material  temperatures  at  the  outer  boundary  of  insulation 
layer,  Tw 4,  and  the  accommodation  factor  of  the  material  tem¬ 
perature,  AF t  =  T(2D)/r(lD),  are  shown  in  Fig.  15a  and  b, 
respectively.  In  general,  AFj  exhibits  a  synergy  with  AF^, 
but  varies  in  the  range  between  0.75  and  0.74  in  the  present 
investigation. 

6.  Conclusions 

This  work  presents  efforts  to  establish  integrated  electro¬ 
chemical  and  thermal  analyses  of  a  planar  SOFC.  The  cell 
performance  is  evaluated  by  combining  an  electrochemical 
model  with  a  simplified  thermal  approach,  which  can  be  im¬ 
plemented  for  various  flow  configurations.  Based  on  ‘quasi¬ 
conduction’  procedures,  a  finite  element  model  is  employed 
to  simulate  the  2D  heat- transfer  behaviour  on  an  interconnect 
with  insulation  envelopes.  This  calculates  not  only  the  tem¬ 
perature  distribution,  but  also  the  heat  flux  dissipated  to  the 


environment.  The  latter  then  serves  as  a  correction  term  to 
the  ECM  for  refined  evaluation  of  cell  performance,  so  that 
a  unified  temperature  distribution  with  a  heat  loss  effect  can 
be  obtained  via  iteration  processes. 

It  has  been  verified  that  a  temperature  accuracy  of  better 
than  0.1%  can  be  achieved  within  a  few  iterations.  In  addi¬ 
tion,  this  work  relieves  the  constraint  of  thermal  boundaries 
for  the  active  area  of  SOFC,  and  realistically  mimics  the  heat- 
transfer  phenomenon  at  the  periphery.  By  combining  ECM 
with  FEM,  cost-intensive  and  time-consuming  CFD  efforts 
are  circumvented,  and  consistent  and  synergetic  thermal  char¬ 
acteristics  are  generated.  Hence,  it  can  be  concluded  that  an 
efficient  computational  methodology  for  integrated  electro¬ 
chemical  and  heat-transfer  analyses  has  been  accomplished, 
which  will  be  utilized  in  the  future  to  acquire  the  necessary 
design  data  for  prototyping  SOFC  system  applications. 

A  bulk  heat-transfer  model  for  simplified  design  analysis 
has  also  been  developed.  By  implementing  various  empiri¬ 
cal  correlations  from  textbooks,  the  model  can  be  utilized  to 
study  either  normal  or  lateral  heat  transfer  phenomena  in  a 
planar  cell,  of  which  the  latter  has  been  the  major  concern  in 
this  study.  In  engineering  practice,  insulation  materials  would 
be  deployed  as  an  enclosure  for  an  SOFC  stack  to  reduce  heat 
loss  to  the  environment.  The  concept  of  total  heat  resistance 
would  then  be  employed  to  facilitate  ID  heat- transfer  analy¬ 
sis,  and  the  relative  magnitude  of  the  resistance  due  to  differ¬ 
ent  heat-transfer  mechanisms  or  material  properties  could  be 
evaluated  accordingly.  Hence,  predominant  parameters  that 
affect  heat-transfer  behaviour  could  be  readily  determined. 
Moreover,  some  ‘accommodation  factors’  to  correlate  the  ID 
results  with  the  counterparts  of  2D  FEA  have  been  deduced 
from  the  presently  available  data.  It  is  expected  that  more 
comprehensive  correlations  can  be  generated  with  extended 
parametric  studies,  which  will  be  beneficial  for  providing 
rapid  prototyping  processes. 
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